Empirical Calculations of the Formation Energies of Intrinsic Defects in Lithium Niobate 
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Abstract 

Results of our shell-model calculations concerning the determination of 
empirical parameters of inter-ionic potentials and the formation ener- 
gies of point defects in LiNb03 are presented. We employed the relaxed 
fitting method which is particularly appropriate for the modeling of low- 
symmetry structure. Attention is paid to guarantee the simulation pre- 
cision of the parameters who are sensitive for defect formation energy 
calculations. It is shown that Li Frenkel disorders are the main intrinsic 
defects in stoichiometric lithium niobate and Li vacancy model is the 
dominant defect species in congruent lithium niobate. 

1 Introduction 

Lithium Niobate (LiNbC>3) is widely used ferroelectric material because of its favorable opti- 
cal, piezoelectric, electro-optic, elastic, photoelectric and photorefractive properties. LiNbC>3 
has a tendency to non-stoichiometry with R=[Li]/[Nb]< 1. Such crystals therefore have very 
high concentration of intrinsic defects. 

The intrinsic defect structures of LiNbC"3 are widely investigated because many physical 
properties of the crystal depend to a certain extent on it. The required local charge neutrality 
can theoretically be guaranteed by one of the models listed in Table 1. 

Table 1. Notations of the intrinsic defect structure of LiNbC>3 



Model No. 


Parametric formula 


Kroger- Vink notation 


Authors 


1 


Li^xNbOg-x 


2V Li V 


Fay|l| 


2 


[Li!_5xNb x ]Nb0 3 


Nb Li 4V Li 


Lernerf^] 


3 


[Li 1 _ 5x Nb 5x ]Nb 1 _4x0 3 


5Nb Li 4V Nb 


Abrahams jD 



Model 1 was suggested by Fay et al. (in 1968) that Li20 deficiency in LiNbOs is accom- 
modated by the formation of Li and O vacancies E]]. Model 2 was suggested by Lerner et al. 



(in 1968) that, in crystals with R= [Li] / [Nb] < 1, excess Nb 5+ ions are substituted for Li + in 
their sites. To satisfy the electro-neutrality condition, there will then arise lithium vacancies. 
The Nb sublattice is full occupied and every four vacant Li site (charge of -4) are statistically 
accompanied by one Li site occupied by a Nb ion (difference charge of +4 at this site)[Q]. 
Model 3 was suggested by Abrahams and Marsh (in 1986), that the O sublattice remains 
100%, where the excessive charge of five Nb +5 ions on Li-site (charge of +20) is compensated 
by four vacant Nb-site (charge of -20) Q . 

Although careful experimental investigation of the intrinsic defect structure have brought 
some insight into a microscopic understanding, there is no general agreement upon any one 
of the defect model 

Based on the shell model ||, Donnerberg et al. carried out pioneer atomic simulation 
studies on defects in LiNbOs crystals || Qj/J, by which they can then compare these defect 
structure models on the same energetic analysis in terms of their formation energies. 

In such empirical approach the model parameters as well as the potential parameters are 
determined in terms of the empirical parameterization, i.e. all the parameters involved are 
adjusted to reproduce as closely as possible the measured crystal data. 

Nevertheless, because of the low symmetry and the covalency cohesion involved for the 
ternary oxide lithium niobate, the practical limitations of the empirical parameterization 
bring about that either the residual strain will be present, if the crystal properties (cohesion, 
elastic constants, dielectric constants and phonon spectrum) are fitted exactly, or conversely 
the poor crystal properties may be obtained if all ion and bulk strain is eliminated. 

Since the fitting result is not unique, an optimum parameter set has to be chosen deliber- 
ately according to succeeding researches. Therefore, we refine upon some parameters which 
are sensitive for defect formation energy calculations, whereas loosen the errors for the other. 
Moreover, the fitting result is concerned with the path of the optimization. In the present 
study, the relaxed fitting approach by GaleQ] is adopted. 

Based on the optimum parameter set, the defect formation energies can then be calculated 
by using the division strategy of Lidiard and Norgett [9] |T(J] incorporating with Mott-Littleton 
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approximation |Tl| . 

Lithium deficient non-stoichiometric LiNbOs can be imaged as the result of Li20 out- 
diffusion from stoichiometric one. There are different out-diffusion modes assumed in respect 
to corresponding defects emerging in the host LiNb03 after the removal of the L^O. By 
comparing three diffusion modes in terms of the energy needed for removing one Li20, one 
can conclude which intrinsic defect model will be energetically favorable in lithium deficient 
non-stoichiometric LiNb03. 



2 Parameterization 



According to the shell model, each ion is assumed to consist of an core (charge X, mass— > 
Mi on ) and an electron shell (charge Y, mass— >0). Therefore the number of constituent units 
of the crystal is doubled. The polarization of the crystal is caused by the displacements of 
ions and the relative displacements between cores and shells. 
There are four kinds of interaction for the crystal. 

(1) The interaction between shell and core belong to the same ion, 

c = ~^k sc \r s — r c | (1) 

in which the spring constant k sc characterizes the harmonic coupling between the shell and 
core. 

(2) The Coulomb interactions between ions i and j, 

XjXj XoY^ YjXj YiYi , , . 

(t>couiomb = | — + i — r^+i — r^ + i — -^—\ { + i ( 2 ) 

\^ic ^jc\ |*"ic ^js\ l^is ^jc\ \^is ^js\ 

where i ^ j emphasizes that there is no Coulomb interaction between shell and core belong 
to the same ions. 

(3) The short-range interactions between shells are described with Buckingham potential 
form: 

(t>non-Couiomb{r) = Mj exp(r / pij) - Qj /r 6 (3) 



where r is the distance between the centers of electron shells of species i and j. Because the 
radius of cation is relatively small the interaction between cations is ignored. Besides, the 
van der Waals terms — Cy/r in anion-cation interaction is omitted too. 

(4) Three body bond-bending terms are used to model the directional character of covalent 
bond: 

<f>bend{0) = \k B (e-d f (4) 

where ks is the bond-bending force constant, is the bond angle of O-Nb-O, and 9$ the 
equilibrium bond angle. 

Accordingly to the equations listed above, there are 14 parameters involved. Since the 
polarizability of Li ion, au = Y^/k^ , is very small, so we take approximately fej£ — > oo, 
Yhi = . Therefore the member of undetermined parameters reduced to 12. A parameter set 
can be regarded as a point in the " 12-dimensional parameter space" . 

We can properly choose an initial parameter set as a starting point and then move the 
point closed to the reasonable position by an optimization procedure. In the present study, 
"relaxed fitting" approach by Gale is employed. This approach based on displacements of 
the energy-minimized configuration from the experimental structure is found to be superior 
to conventional fitting in which zero gradients of the experimental geometry is the aim. The 
relaxed fitting yields the exact displacements and genuine physical properties and appears to 
be particularly appropriate for the modeling of low-symmetry structure ||. 

In order to quantify the degree of the optimization, we define an objective function as 
following: 

F = U^2(S iCa lc - SiExpr) 2 / Sf Ref + V^(P jC alc ~ PjExprf 7 P jRef u i v > ( 5 ) 
i j 

It is obviously that F > and small F implies small fitting errors or high degree opti- 
mization. {S} and {P} denote the quantities of crystal structure and property respectively. 
u and v are weighting factors. Generally, the precision of reproduced structure should be 
more emphasized in comparing with that of reproduced properties. Therefore, in the present 
study, we take u = 100 and v = 1. 
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(SiCaic ~ SiExpr) 2 / Sf Ref and (PjCaic - Pj Ex P r) 2 / P jR e f represent the fitting errors between 
that calculations and experiment date. In order to avoid the irrationality resulting from small 
denominator, we take the referential values instead of the experiment data in some cases. For 
instance, the fractional coordinate y(=0.3446) of the oxygen ion in hexagonal cell is many 
times larger than x(=0.0492) and z(=0.0647). But in practice, we take their referential 
values as ^-R e f=jRef =z Ref = ^-^- Similar treatments are adopted for elastic constants(we use 
10xl0 10 N m -2 as their referential values). 

The object function, F, is then minimized in terms of Davidon-Fletcher-Powell variable- 
metric algorithm. Since the optimization often gets stuck in local minimum, and so we have 
to try hundreds of initial parameter sets and then choose the globe minimum from hundreds 
of local minima. 



Table 2. The two parameter sets of ionic interaction and shell model 



Parameters of short- 


Set I 


Set II 


Parameters of shell 


Set I 


Set II 


range interaction 






model 






^Nb-o/ev 


5121.9 


6373.6 


*W|e| 


-36.0 


-36.0 


PNb-o/A 


0.27543 


0.27211 


C7(evA- 2 ) 


16553. 


13374. 


Au-o/ev 


144.27 


246.64 


Yo/\e\ 


-2.5376 


-2.4660 


PLi-O /A 


0.35614 


0.32554 


O(evA-) 


11.397 


11.459 


Ao-o/ev 


270.35 


198.80 


Parameters of bond- 
bending 


Set I 


Set II 


Po-o/A 


0.40510 


0.42612 


fcs/(ev rad~ 2 ) 


0.0 


0.51548 


Co-o/evA«) 


11.117 


21.277 


v° 




90 



Table 2 lists two parameter sets obtained, their corresponding fitting results upon crystal 
structure and crystal property are listed in Table 3 and 4 respectively. 



Table 3: Comparison between experimental and calculated values of the structure of 

LiNbOs crystal in hexagonal cell. 
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Structure 


Exptl.|lJ 


Calcd. by set I 


Calcd. by set II 


a 


5.148 


5.1477 


5.2547 


c 


13.863 


13.564 


13.906 


XOc 


0.0492 


0.01481 


0.01816 


XOs 


0.0492 


0.06293 


0.06713 


YOc 


0.3446 


0.37173 


0.35736 


yos 


0.3446 


0.34061 


0.36516 


ZOc 


0.0647 


0.06385 


0.07553 


ZOs 


0.0647 


0.07164 


0.08350 


z Nbs 


0.0 


-0.00026 


0.00100 


ZLi 


0.2829 


0.28738 


0.27500 



where a and c are lattice constants, (x, y, z) are fractional coordinates. The subscript 
c and s represent ion core and electron shell respectively. We consider the core and 
the shell are not coincide when crystal is under equilibrium status, i.e. there is a net 
dipolc within the unit cell. The core and shell of Nb ion are in the triadaxis with 
coordinates (0, 0, 0) and (0, 0, ZNbs) respectively. Because the polarization of Li ion 
is omitted its core and shell has the same coordinate (0, 0, Zli). 

Table 4: Comparison of experimental and calculated data of crystal properties for LiNb03 



Crystal 
properties 


Exptl.|13] 


Calcd. of 
Tomlinson 


Calcd. of 
our set I 


Calcd. of 
our set II 


Cn 


20.3 


25.7 


21.30 


21.27 


Cl2 


5.3 


12.5 


8.91 


5.14 


Cl3 


7.5 


10.6 


6.90 


7.47 


Cl4 


0.9 


-2.2 


-3.67 


2.66 


C33 


24.5 


27.6 


22.28 


25.19 


C44 


6.0 


8.8 


6.92 


6.86 


C66 


7.5 


6.6 


6.20 


8.06 


£11 (0) 


84.1 


29.7 


81.07 


27.48 


£33(0) 


28.1 


61.8 


33.56 


40.74 


£ll(oo) 


46.5 


2.06* 


30.45 


17.68 


£33(00) 


27.3 


1.83* 


21.35 


25.68 



where £n(0), £33(0) are low frequency dielectric constants, £11(00), £33(00) are high 
frequency dielectric constants, Cn, C12, C13, C14, C33, C44, C66 (W 10 Nm~ 2 ) elastic 
constants. 

It seems that the fitting results of the crystal properties result from set I are better than 
that of set II. However the parameters set I leads to imaginary modes of phonons, which 
indicates structural unstable. Therefore parameters set I is given up and parameter set II is 
the accepted result to be used for intrinsic defect studies. 

In comparing with the results of Tomlinson et aZ.|| (see Table 4), two of the "most 
noticeable deviations", as indicated by Tomlinson et al., namely the static dielectric constant 
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£33(0) and the elastic constant C14 are reduced in the present study. Moreover, according 
to the assessment of us, the calculated high frequency dielectric constant en (00) and e33(oo) 
results from the parameter set reported by Tomlinson et al. would be 2.06 and 1.83(marked 
with '*' in table 4), such deviations with magnitudes in order, are diminished considerablely 
as well in our study. 



3 Lattice formation energy and point defect formation energy 



The lattice formation energy of LiNb03 is deduced from parameter set II listed in Table 5. 
Then we extrapolate parameter set II to the empirical calculation for the lattice formation 
energies of Li 2 and Nb 2 Os as well. The results obtained are also listed in Table 5. 

By using the division strategy of Lidiard and Norgett incorporating with Mott-Littleton 
approximation, we can then calculated the point defects in LiNb03. The formation energies 
of different point defects listed in Table 5 are calculated by using an improved homemade 



program similar to HADES [ 10 The results are then combined to give Prankel and 
Schottky energies which are also represented as the total energy and the energy per defect 
(shown in parenthesis, i.e. the Frenkel-pair energy divided by two, the energy of the Schottky 
quintet in LiNb03 divided by five, the energy of the Schottky trio in Li 2 divided by three 
and the energy of the Schottky heptad in Nb 2 Os divided by seven), and are shown in Table 
6. 



Table 5: Some typical lattice energies and defect energies 



Lattice 
energies 


LiNb0 3 


E(LiNb0 3 ) Mni4 ce u 


-183.95 ev 


Li 2 


E(Li 2 0) unit ce ii 


-31.67 ev 


Nb 2 5 


E(Nb 2 05) nn j£ ce u 


-328.10 ev 


Formation 
energies 
of defect 


Li vacancy 


ECVJa) 


8.73 ev 


Li interstitial site 


E(LiJ) 


-6.76 ev 


Nb vacancy 


E(V^ b ) 


129.46 ev 


Nb interstitial site 


E(Nbf) 


-119.15 ev 


vacancy 


E(V£) 


18.83 ev 


interstitial site 


E(0'/) 


-10.16 ev 


Nb antisite 


E(Nbg) 


-114.82 ev 



7 



Table 6 The formation energies of Frenkel and Schottky defects in L1NDO3 crystal. 



Li Frenkel pair 


1.97 (0.98) ev 


Nb Frenkel pair 


10.31 (5.15) ev 


Frenkel pair 


8.67 (4.33) ev 


LiNb0 3 Schottky 


10.73 (2.15) ev 


Li 2 Schottky 


4.62 (1.54) ev 


Nb 2 5 Schottky 


24.97 (3.57) ev 



4 Intrinsic defect structure in Li deficiency LiNb03 



Li 2 deficiency non-stoichiometric LiNb03 could be imaged as the result of Li 2 out-diffusion 
from stoichiometric LiNbO,3. Theoretically there are three Li 2 out-diffusion modes related 
to three intrinsic defect models respectively, namely: 

'LiNb0 3 « — ► Li 2 + 2Vl; + V" (model 1) (6) 



'LiNb0 3 « — > 3Li 2 + 4V' Li + Nb£* (model 2) (7) 



'LiNb0 3 * — ► 3Li 2 + 4V^j' b + 5Nb£- (model 3) (8) 

The niobium vacancy (V^j b ) with effective charge of — 5|e| tends towards to the niobium 
antisite (Nb^*) with effective charge of +4|e|, which decreases their total formation energy. 
Table 7 illustrates the formation energy dependence of the distance between V^ b and Nb^*. 
It can be seen that the nearest neighboring V^-Nb^* pair, complex defect [V^Nd^*], has 
the lowest formation energy. Therefore, a slight modification of (8) yields 

' LiN b0 3 <— > 3Li 2 + 4(vg b Nb&) / + Nb^' (9) 



Table 7: The relation between V^j b — Nb^* distance and formation energy. 



V^ b — Nb^* distance 


infinite 


next-next 
nearest 
neighbor 


next 
nearest 
neighbor 


nearest 
neighbor 


total formation energy 


14.64ev 


12.73ev 


12.29ev 


9.98ev 



8 



Table 8 contains our calculation results in which all the energies needed for (6)~(9) out- 
diffusion modes are converted to per L^O. In order to make a comparison with the previous 
studies, the corresponding results reported by Donnerberg et al.[|j are listed in Table 8 as 
well. 



Table 8 The defect energies per Li20 of four typical defect model. 



Process of out-diffusion 


Eq. (6) 


Eq. (7) 


Eq. (8) 


Eq. (9) 


Calcd. by us 


4.62ev 


3.01ev 


10.89ev 


4.68ev 


Calcd. by Donnerberg 


5.82ev 


4.56ev 


15.2 ev 


10.7ev 



The energy difference between Eq.(6) and Eq.(7) indicates the model 1 (oxygen vacancy 
model) is energetically unfavorable. Though the energy denoted by Eq.(9) is much less than 
corresponding result reported by Donnerberg et al. [4], it is still 1.5 times larger than the 
energy denoted by Eq.(7). 

The fact that the energy denoted by Eq.(7) is the lowest implys the defect structure model 
2 (lithium vacancy model) is energetically favorable. 



5 Conclusion and discussion 



The simple defect model 1 was rejected by several experimental evidences. Until 1993 the 



model 3 was favored. However, the earlier NMR studies[15] which had strongly supported 



the niobium vacancy model is in doubt[16|. 

The x-ray single-crystal diffraction and the TOF neutron powder diffraction investigation 
by Iyi et al. strongly supported the Li-site vacancy model]!?]]. Wilkinson et al. examined the 
congruently melting lithium niobate samples by powder x-ray diffraction techniques and found 



the Nb site is fully occupied [18]. The x-ray and neutron powder diffraction experiments by 
Zotov et al. shown that the Li-sites vacancy model describes best the average defect structure 



in congruent lithium niobate[19|. 

Recently, by using a formula for calculating the soft mode frequency and Curie temper- 
ature to the analysis of the defect structure, Safaryan et al. concluded that the lithium 
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vacancy model best describes the defect structure of nonstoichiometric lithium niobatcpC|1. 

According to the present empirical study, we can conclude that lithium deficiency of the 
LiNb03 is accompanied mainly with the formation of lithium vacancies and niobium antisite. 

In comparing with the result reported by Tomlinson et al., the degree of fit between the 
calculated crystal properties and the observed data has been improved in the present study, 
especially the high frequency dielectric constants. Even though there exist still obviously the 
disagreement between En(oo) calculated and experimental values, it will take less effect for 
the polarization behavior of the ion under the electric displacement (D), which is important 
for successive studies of intrinsic defects. The reason is that the polarization density P can 
be expressed as P = ^j^D (deduced from D = sqE + P = sqeE), and so the proportionality 
factors are 0.944 and 0.978 corresponding to en(oo) calculated value of 17.1 and experimental 
value of 46.5 respectively. 
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